rm(list=ls())

load("TSCS_data.RData")

# Figure E.5

models_plm <- list()
world_regions <- levels(as.factor(panel$region))

for (world_regions in 1:17) {
  panel %>% filter(region!= world_regions)  %>% 
    as.data.frame() %>% 
    pdata.frame(index=c("country_name","year"),row.names = FALSE) ->
    data
  
  model <- plm(
    citizen_support ~ lag(citizen_support,1:2)  +
      antipluralism_gov_seat_share_lag + 
      lag(polarization,1) +
      lag(clientelism_gov_seat_share,1) + 
      lag(region_citizen_support,1) + 
      lag(region_antipluralism_gov_seat_share,1) + 
      lag(region_liberal_democracy,1) + 
      lag(democratic_stock,1) + 
      lag(GDP_capita,1) + 
      lag(life_expectancy,1),
    data = data,
    model = "within",
    effect = "twoways")
  
  models_plm[[world_regions]] <- model
}

coeffs_lm <- vector()
se_lm <- vector()

for (i in 1:17) {
  coeffs_lm[i] <- models_plm[[i]][["coefficients"]][["antipluralism_gov_seat_share_lag"]]
}

vcm_lm <- lapply(models_plm, function(x) plm::vcovBK(x, cluster=c("group")))
rse_lm <- lapply(vcm_lm, function(x) sqrt(diag(x)))

for (i in 1:17) {
  se_lm[i] <- rse_lm[[i]][["antipluralism_gov_seat_share_lag"]]
}

coeff_df <- as.data.frame(cbind(coeffs_lm,se_lm))
coeff_df %>% mutate(margin_of_error = se_lm*1.65,
                    CI_lower = coeffs_lm - margin_of_error,
                    CI_upper = coeffs_lm + margin_of_error) %>%
  mutate(model = row_number())->
  coeff_df

pdf(file = "fig_e6a.pdf",width = 5,height = 9)

ggplot(coeff_df,aes(x=coeffs_lm,y=model))  + 
  geom_errorbarh(aes(xmin=CI_lower, xmax=CI_upper), colour="grey") +
  theme() + xlab("Coefficient") + 
  ylab("Model")  +
  geom_vline(xintercept = 0) + 
  geom_point() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 

dev.off()

# GMM

models_pgmm <- list()

for (world_regions in 1:17) {
  panel %>% filter(region!= world_regions)  %>% 
    as.data.frame() %>% 
    pdata.frame(index=c("country_name","year"),row.names = FALSE) ->
    data
  
  model <- pgmm(
    citizen_support ~ lag(citizen_support,1:2)  +
      antipluralism_gov_seat_share_lag + 
      lag(polarization,1) +
      lag(clientelism_gov_seat_share,1) + 
      lag(region_citizen_support,1) + 
      lag(region_antipluralism_gov_seat_share,1) + 
      lag(region_liberal_democracy,1) + 
      lag(democratic_stock,1) + 
      lag(GDP_capita,1) + 
      lag(life_expectancy,1)|
      lag(downturn, 3:5),
    data,
    effect = "individual",
    model = "onestep",
    transformation = "ld",
    indexes = c("country_name", "year"),
    robust = T)
  
  models_pgmm[[world_regions]] <- model
}

coeffs_gmm <- vector()
se_gmm <- vector()

for (i in 1:17) {
  coeffs_gmm[i] <- models_pgmm[[i]][["coefficients"]][["antipluralism_gov_seat_share_lag"]]
}

vcm_gmm <- lapply(models_pgmm, function(x) plm::vcovHC(x, cluster=c("group")))
rse_gmm <- lapply(vcm_gmm, function(x) sqrt(diag(x)))

for (i in 1:17) {
  se_gmm[i] <- rse_gmm[[i]][["antipluralism_gov_seat_share_lag"]]
}

coeff_df <- as.data.frame(cbind(coeffs_gmm,se_gmm))
coeff_df %>% mutate(margin_of_error = se_gmm*1.65,
                    CI_lower = coeffs_gmm - margin_of_error,
                    CI_upper = coeffs_gmm + margin_of_error) %>%
  mutate(model = row_number())->
  coeff_df

pdf(file = "fig_e6b.pdf",width = 5,height = 9)

ggplot(coeff_df,aes(x=coeffs_gmm,y=model))  + 
  geom_errorbarh(aes(xmin=CI_lower, xmax=CI_upper), colour="grey") +
  xlab("Coefficient") + ylab("Model")  +
  geom_vline(xintercept = 0) + geom_point() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 

dev.off()
